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POLAR  2.0  VALIDATION  AND  PREFLIGHT  SPEAR  I  CALCULATIONS 


G.  A.  Jongeward.  M.  J.  Mandell,  J.  Lilley.  and  I.  Katz 
S-CUBED  Division  of  Maxwell  Laboratories,  inc. 


1.  INTRODUCTION 

In  this  report,  we  present  a  series  of  calculations  used  to  test  and  validate  the 
POLAR  2.0  code.  The  goal  in  this  first  limited  validation  phase  of  POLAR  2.0  is  to 
model  the  electron  collection  experiment.  SPEAR  I  (Space  Power  Experiments  Aboard 
Rockets).  The  POLAR  code  (hereinafter  referred  to  as  POLAR  1.1)  was  developed  by 
S-CUBED  for  the  Air  Force  Geophysics  Laboratory.  It  is  currently  in  the  second  year 
of  a  3-year  validation  phase.  Concurrent  with  this  validation  effort,  the  POLAR  2.0 
code  is  being  developed.  POLAR  2.0  contains  the  same  physical  models  as 
POLAR  1.1  but  allows  multiple  grids.  This  increased  versatility  is  required  when  both 
a  large  computational  grid  and  localized  high  resolution  are  needed.  Examples  of 
these  situations  are:  SPEAR  I.  which  has  important  physical  details  on  the  centimeter 
length  scale  and  produces  plasma  disturbances  meters  from  the  rocket;  and  the 
charging  of  an  astronaut  in  the  wake  of  a  polar  orbiting  shuttle  in  which  the  details 
of  the  astronaut's  suit,  extra-vehicular  maneuvering  unit,  and  the  high  voltage  shuttle 
charging  must  all  be  resolved. 

The  report  is  organized  as  follows:  Section  2.  contains  an  overview  of  POLAR  2.0 
and  a  brief  description  of  the  data  structure  with  a  comparison  to  POLAR  1.1.  In 
Section  3..  we  give  a  review  of  the  physical  models  in  POLAR  2.0  that  are  relevant 
to  the  calculations  presented  here.  Section  4.  contains  test  calculations  for  a  10-cm 
sphere  in  a  magnetic  field.  Space  charge  effects  on  potentials,  magnetic  field  effects 
on  particle  trajectories,  and  self-consistent  magnetosheath  calculations  are  presented. 
Section  5.  contains  a  description  of  the  SPEAR  I  object,  calculational  grid  structure, 
and  the  results  of  two  preflight  calculations  of  current  collection  by  the  SPEAR  I 
experiment.  Section  6.  compares  POLAR  2  0  and  NASCAP/LEO  calculations  of  the 
SPEAR  I  current  collection. 
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2.  OVERVIEW  OF  POLAR  2.0 

2.1  Introduction 

The  original  POLAR  code.  POLAR  1.1.  has  been  extended  to  include  nested  grids. 
Nested  grids  allow  calculations  at  higher  surface  potentials  without  loss  of  object 
surface  resolution.  Additionally,  calculation  size  limitations  have  been  removed.  The 
approach  taken  was  to  modify  the  computational  algorithms  of  POLAR  1.1  as  little  as 
possible  by  adding  multiple  grid  capabilities  at  the  data  handling  level.  Changes  were 
made  in  the  general  data  handling  routines,  the  file  handlers,  and  in  the  calculation 
modules  to  deal  with  grid  interfaces. 

2.2  General 

General  data  handling  has  been  modified  in  two  areas:  the  method  of  storing  data  on 
the  MSIO  files,  and  the  use  of  the  PUFF  array.  The  new  MSIO  routines  can  be 
used  by  POLAR  1.1.  POLAR  2.0  uses  multiple-indexed  MSIO  files.  When  OPENMS 
is  called  to  read  the  file  index,  the  number  of  keys  to  be  used  to  address  the  file 
records  is  defined.  File  19.  the  data  file  containing  quantities  referred  to  by  name. 
SRFV  and  so  on.  adds  a  second  key  for  the  grid  number. 

The  MRBUF  routine,  which  in  POLAR  1.1  sets  up  the  VPROPS  common  block,  has 
been  moved  to  the  PATCHV  module.  The  location  of  POLAR  data  in  both  the 
MSIO  data  files  and  the  PUFF  array  is  determined  during  object  definition.  This  has 
removed  many  of  the  limits  on  the  grid  size.  The  discussion  of  the  PATCHV  module 
describes  these  modifications  in  more  detail. 

The  complexity  of  the  data  structure  used  in  POLAR  2.0  precludes  the  use  of 
statement  functions  to  calculate  data  locations.  The  statement  functions.  IPUFAD 
and  IPUFZS.  have  been  replaced  by  integer  functions  of  the  same  names.  Two  new 
data  vectors.  SGRD  and  NSGD.  have  been  added  to  describe  the  elements  and  nodes, 
respectively,  of  each  grid. 

2.3  Object  Definition 

The  main  difference  between  the  two  versions  of  POLAR  is  apparent  to  the  user 
during  object  definition.  Because  POLAR  2.0  is  oriented  towards  multiple  grid 
objects,  the  definition  of  the  grids  and  the  objects  within  them  is  much  more 
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complicated  than  in  POLAR  1.1.  The  grid  structure  presently  in  POLAR  2.0 
presumes  the  eventual  addition  of  generalized  geometries  to  POLAR.  Grids  may  be 
empty  or  contain  ORIENTed  VEHICL  denifitions.  The  outermost  grid  is  used  to 
define  the  mesh  size  in  meters  for  all  of  the  inner  grids  and  is  also  the  only  one  that 
may  be  stepped  as  the  POLAR  1.1  calculational  grid. 

Grids  are  subject  to  the  following  restrictions: 

1.  Grids  must  be  defined  in  order  of  nesting,  outer  first. 

2.  Grids  must  be  completely  enclosed  by  any  outer  grids,  although  boundaries 

may  be  shared. 

3.  The  grid  origin  must  lie  on  a  node  in  the  outer  grid. 

4.  The  ratio  of  the  inner  to  outer  mesh  size  must  be  an  integer  greater  than 

zero. 

5.  Grid  boundaries  must  be  on  nodal  planes  in  the  outer  grid  in  order  to 
interpolate  boundaries  values.  This  means  the  grid  length  must  be 
evenly  divisible  by  the  mesh  ratio  of  the  grid  to  the  outer  grid. 

6.  Inner  grids  may  not  touch  other,  equal  level  inner  grids.  (Not  checked 
for.) 

Current  restrictions  to  prevent  surface  subdivision  are  as  follows: 

1.  Grid  boundaries  cannot  touch  nonzero  element  table  entries  in  the  outer 
grids. 

2.  A  grid  cannot  contain  objects  defined  in  an  outer  grid. 

In  order  to  help  the  user  maneuver  objects  into  their  desired  locations  in  the  mesh, 
several  enhancements  have  been  made  to  ORIENT.  Grids  may  be  expanded  and/or 
shifted. 

A  new  module  was  necessary  for  mesh  and  object  definition  for  POLAR  2.0 
calculations.  This  module.  PATCHV.  combines  any  number  of  object  definitions 
created  by  the  POLAR  1.1  modules.  VEHICL  and  ORIENT,  with  empty  grids  to 
produce  the  final  calculation  mesh.  Presently.  PATCHV  is  an  interactive  module  that 
builds  the  mesh  from  the  outer  grid  into  the  inner  grids.  Since  all  of  the  grid 
information,  such  as  nx.  ny.  and  nz,  and  the  number  of  surfaces  are  known  during 
mesh  definition,  the  MSIO  Files  11.  9,  and  19  are  dynamically  defined.  This  means, 
for  example,  that  the  POT  z  slices  will  probably  be  found  at  different  index  key 
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locations  for  different  mesh  definitions.  A  permanent  stack  structure  is  defined,  which 
contains  mesh  and  vector  property  information.  The  locations  of  pointers  and  vector 
definition  will  also  vary  from  run  to  run. 


Fig.  2.1  shows  the  pointer  and  data  structure  for  POLAR  2.0.  The  figure  shows  the 
information  as  it  is  placed  in  the  common  block  CBUF.  The  information  in  CBUF 
breaks  into  four  groups.  These  are:  1.  the  pregrid  information:  2.  the  pointers  to  the 
grid  properties:  3.  the  vector  properties;  4.  the  data. 
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The  first  group  is  the  pregrid  information.  This  region  contains  global  information 
such  as  number  of  grids  and  mesh  size.  It  also  contains  pointer  to  the  start  of 
information  on  the  grids  in  Region  2  as  shown  in  Fig.  2.1.  As  indicated.  Region  2 
starts  at  location  igridO+1  and  contains  both  values  and  pointers  to  the  grid 

properties.  The  information  is  grouped  by  grid,  and  the  grid  data  follows  immediately 

each  grid's  pointer  information.  Region  3  contains  the  vector  property  information. 

The  starting  point  for  this  information  is  iavecO+1.  The  vector  properties  are  stored 
by  grid,  so  the  total  length  of  this  region  is  (number  of  grids)  x  (number  of  vectors) 
x  (number  of  vector  properties).  This  region  contains  the  IAVEC  data  of  POLAR  1.1. 
The  majority  of  the  information  on  the  vectors  is  actual  values.  A  few  pointers  are 

mixed  in  with  this  information:  for  example,  the  starting  location  of  the  vector  data  in 

PUFF  is  recorded  in  the  second  vector  property.  The  final  region  is  for  the  paged 
information.  This  is  where  all  the  information  read  and  written  to  disk  by  POLAR  s 
data  base  manager  is  kept  in  memory. 

2.4  Calculations 

The  POLAR  1.1  NTERAK  module  has  been  divided  into  two  separate  executables  in 
POLAR  2.0  for  development  purposes:  POT  performs  the  PWASON  and  charge 
calculations,  while  DENS  calculates  initial  geometric  ion  densities  and  the  "DIONS" 
found  by  "CURREN” 

At  the  present  level  of  development.  POLAR  2.0  is  not  really  a  very  user-friendly 
code  It  has  been  tested  with  respect  to  the  algorithms  used  to  handle  data  sharing 
between  grid  boundaries,  loops  over  the  grids,  and  general  soundness  of  the  data 
handling.  Work  needs  to  be  done  to  speed  up  calculations.  Several  physical  features 
found  in  POLAR  1.1  have  not  been  fully  implemented  yet 

3.  POLAR  2.0  ELECTRON  MAGNETOSHEATH 

In  this  section,  we  review  the  physical  models  and  calculational  methods  that  go  into 
a  self-consistent  calculation  of  a  magnetosheath.  We  will  first  discuss  a  semianalytic 
approach  to  calculate  the  potentials  that  is  valid  for  symmetric  objects.  This  will  be 
followed  by  a  description  of  the  electron  particle-pushing  algorithms  and  finally  the 
computation  of  a  self-consistent  magnetosheath. 
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3.1  Symmetric  Potentials 

The  potential  solver  in  POLAR  2.0  solves  Poisson's  equation. 


The  charge  density  can  result  from  a  geometrical  ion  calculation,  barometric  law, 
analytical  formulae,  or  from  trajectory  tracking.  For  high  voltage  electron  attracting 
sheaths,  wakes  do  not  produce  a  significant  change  in  the  sheath.  As  a  starting 
point  for  a  calculation.  POLAR  2.0  uses  an  analytic  formula  relating  the  charge 
density  to  the  potential,  p  =  p(f. E)  This  analytic  relationship  produces  the  correct 
sheath  size  and  shape  for  symmetric  geometries  (spherical,  cylindrical,  planar)  and 
gives  a  reasonable  starting  point  for  more  complex  shapes. 


The  analytical  form  for  p(f. E) 

Rsheath 

for  the  convergence  factor,  - ,  namely, 

r 

,2 


fR 


sheath 


=  const  E*VC 


and  then  choosing  b  and  c  to  best  approximate  the  exact  relationship  in  a  least 
squares  sense.  The  following  form  was  found  to  be  an  excellent  fit  over  a  wide 
range  of  temperatures,  densities,  and  potentials: 


Pit, E)  *  PQ 229 


/ 

a 

UJ 

1.262 

-.509 

r  e 

l  6  . 

.0. 

.Atf. 

.5 


3.2  Electron  Particle  Pushing 

The  Earth's  magnetic  field  strongly  influences  the  motion  of  electrons.  A  thermal 
electron  has  a  gyroradius  of  2.5  cm.  This  is  much  smaller  than  typical  object  and 
sheath  dimensions,  and  much  smaller  than  the  10-m  SPEAR  I  rocket  and  sheath. 
Hence,  electrons  in  regions  of  weak  electric  field  will  ExB  drift  and  in  regions  of 
strong  fields  will  execute  complex  trajectories,  bouncing  many  times  within  the  sheath. 
Some  electrons  will  enter  noncollecting  "trapped”  orbits,  some  will  eventually  be 
collected,  and  still  others  will  leave  the  sheath.  For  these  reasons,  POLAR  11  and 
POLAR  2.0  push  electrons  differently  than  ions.  In  some  instances,  step  pushing  is 
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used:  while  in  other  instances,  a  guiding  center,  "drift"  approximation  is  used. 
Electrons,  like  ions,  are  pushed  using  a  super  particle  that  has  a  current. 

I  =  A  .  .  J  .  ..  For  each  cell,  the  contribution  to  the  plasma  density  by  a 

superparticle  is  <5Q  =  It.  where  r  is  the  time  for  the  particle  to  cross  the  cell.  First, 
we  discuss  the  standard  step  pushing  method,  followed  by  a  discussion  of  the  drift 
method. 

The  equation  of  motion  of  an  electron  in  an  electric  field  is 

m  —  =  eVf  -  eVx§ 
dt 

In  zones  where  the  cyclic  motion  of  the  electron  is  comparable  to  a  zone  size,  step 
pushing  is  used.  A  time  step  that  will  produce  10  steps  to  cross  a  zone  is 
estimated.  The  velocity  and  position  are  updated  using 

$i+1  =  r  ♦  q£ 1  fit 
Xi+1  =  r  ♦  $!fit 

The  rotation  of  velocity  by  the  magnetic  field  is  computed  to  second  order  in  fit  as. 

2  2 

5=5.5  .  [5j  (}*§)„§  §i- 

Using  this  algorithm  for  the  magnetic  field,  rotation-produced  electron  gyroradii 
accurate  to  within  a  percent.  This  2-step  approximation  is  much  quicker  than  the 
more  accurate  matrix  multiple  due  to  the  time  required  to  construct  and  rotate  the 
matrix. 

Resolving  the  drifting  motion  of  electrons  in  the  weak  electric  field  regions  is  very 
CPU  intensive.  To  speed  the  numerics.  POLAR  uses  a  drift  approximation  that 
breaks  the  motion  into  two  components.  The  velocity  and  acceleration  parallel  to  B 
is  unchanged:  while  perpendicular  to  B,  the  velocity  is  the  drift  velocity  with  no 
acceleration  Specifically,  for  the  parallel  velocity. 


m 


dt 
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and  for  the  drift  velocity, 


The  precession  of  V  about  B  is  computed  in  the  drift  frame  where. 


dt 


where  IWI  =  ~~  is  the  cyclotron  frequency. 


This  drift  approximation  is  used  when  the  gyro  motion  is  small  compared  to  the  zone 
size.  When  this  is  the  case,  resolving  the  gyro  motion  would  not  produce  optimal 
accuracy  to  the  trajectory  because  the  local  fields  are  not  precisely  known.  The 
choice  of  when  to  use  the  drift  approximation  is  made  by  comparing  the  local  gyro 
radius  to  the  zone  size  and  switching  to  drift  motion  when  the  gyro  radius  is  one-half 
of  a  mesh  length. 


3.3  Self-Consistent  Sheaths 

The  sequence  of  calculations  used  to  compute  a  self-consistent  electron  magnetosheath 
is  listed  below: 


Step  1 

Compute  the  potentials  approximating  the  space  charge  density  ignoring  magnetic  field 
effects. 


V2*0  =  - 


»(eV) 


Step  2 

Track  electron  trajectories  from  the  sheath  edge  until  one  of  the  following  happens: 
they  hit  the  object;  they  exit  the  sheath;  or  they  enter  a  trapped  orbit.  Use  either 
the  step  push  or  drift  methods  depending  on  the  field  strengths.  Use  the  fields  from 
the  latest  potential  solution. 


^  =  q  t' 

dt 


q  U 


=  y 

dt 


n 
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Step  3 

Using  the  time  required  for  a  particle  to  cross  a  cell.  r.  increment  the  local  plasma 
density  for  use  in  the  next  iteration  of  the  potential  solver. 


<5q 


i 


^sheath^thT 


where  Asheath  is  the  sheath  area  the  superparticle  represents. 


Step  4 

After  tracking  all  superparticles  from  the  sheath  edge,  recompute  the  potentials. 


i  +1 


put'}) 

eo 


where  />({/})  represents  the  plasma  density  functional  of  the  potentials  at  iteration  i 
computed  by  particle  tracking.  We  then  loop  back  to  the  particle  tracking  phase. 
Step  2.  and  iterate  until  the  collected  current  becomes  stable.  For  most  cases, 
stability  can  be  defined  in  terms  of  the  collected  current 


4.  POLAR  2.0  ELECTRON  COLLECTION  VALIDATION  TESTS 


To  validate  the  current  and  potential  algorithms  in  POLAR  2.0.  we  used  as  a 
standard  test  object  a  10-cm  sphere  biased  to  11.000  volts.  The  environment 
consisted  of  a  0.1-eV  plasma  with  density  lO^m  3.  This  was  used  as  a  standard 
test  object  in  the  NASCAP/LEO  code  and  in  a  2  1/2-dimensional  code.  The 
2  1/2-dimensional  code  used  cylindrical  symmetry  to  produce  results  that  are  valid  for 
three  dimensions,  including  magnetic  field  effects.  We  will  present  comparisons 
between  POLAR  2.0  and  these  two  codes. 

The  test  object,  shown  in  Fig  4.1.  is  a  QSPHERE  (quasisphere)  with  a  diameter  of 
three  mesh  units.  It  has  26  faces  composed  of  squares,  rectangles,  and  triangles. 

This  object  is  embedded  in  a  nest  of  three  grids  which  mesh  sizes  0.536  m.  0.134  m. 
and  0.067  m  as  shown  in  Fig.  4.2.  The  outer  grid,  Grid  1,  is  15  by  15  by  15  grid 
units.  Grid  2  is  12  by  12  by  12  grid  2  units,  or  3  by  3  by  3  grid  1  units.  Grid  3. 
which  contains  the  object,  is  10  by  10  by  10  grid  3  units,  or  5  by  5  by  5  grid  2 
units.  The  QSPHERE  is  located  at  (4.  4,  4)  in  grid  3  units.  Three  grids  were  used 
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to  test  the  grid  nesting  in  POLAR  2.0.  The  grid  ratios  were  chosen  to  be  identical 
to  the  grid  ratios  used  in  the  SPEAR  I  calculations. 


\ 

\\ 


Fig.  4.1 
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Fig.  4.2 
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4.1  Space  Charge  Sheaths  without  Magnetic  Field  Effects 

For  this  test,  potentials  were  computed  using  the  analytical  charge  density.  p( E.f). 
during  the  first  iteration  and  trajectory-determined  charge  density  for  subsequent 
iterations.  First,  we  discuss  a  particle  collection  test  which  checked  POLAR's  ability 
to  track  electrons  in  combined  magnetic  and  electric  fields.  This  will  be  followed  by 
a  discussion  of  a  self-consistent  computation  of  the  potentials. 

Fig.  4.3  shows  the  potential  contours  computed  by  POLAR  2.0  using  the  analytic 
charge  density;  therefore,  these  potentials  do  not  include  magnetic  field  effects.  Three 
contours  are  of  particular  interest.  The  first  contour  is  the  conventional  sheath 
contour  0.67*kT.  The  second  is  the  0.75-volt  contour  computed  by  POLAR  as  the 
best  estimate  of  the  sheath  edge.  The  third  is  the  22- volt  contour  that  corresponds 
the  potential,  which  will  drop  to  zero  due  to  spacecharge  shielding  in  one  0.536-m 
zone. 


Fig.  4.3 
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The  analytical  value  for  the  sheath  radius  is  3.6  m.  The  radii  of  the  three  contours 
discussed  above  are  2.9  m,  3.4  m,  4.0  m.  respectively.  The  0.67*kT  contour  is  much 
too  large  and  would  have  been  even  larger  had  the  grid  been  larger.  The  POLAR 
default  0.75-volt  contour  appears  to  be  the  best  value  for  the  sheath  location.  The 
22-volt  contour  is  more  than  one  zone  too  small. 


4.2  Magnetic  Field  Effects  on  Electron  Collection 

Fig.  4.4  shows  trajectories  for  20  electrons  started  3.1  m  from  the  10-cm  sphere. 
The  purpose  of  this  study  was  to  check  the  particle  tracking  of  electrons  in  a 
magnetic  field. 
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As  discussed.  POLAR  uses  both  a  guiding  center  and  a  full  particle  tracking, 
depending  on  the  electric  and  magnetic  field  strengths.  The  magnitude  of  the 
magnetic  field  was  taken  as  0.3  Gauss,  which  gives  a  thermal  electron  a  gyroradius  of 
2.5  cm.  This  is  much  smaller  than  the  mesh  size  of  53.6  cm  in  the  outer  grid. 

Even  for  the  innermost  grid,  this  is  one-third  of  a  zone  Table  4.1  below  shows  the 
input  parameters  for  this  SHONTL  2.0  run.  The  particles  were  started  on  an  arc 
3.1  m  from  the  sphere,  which  are  located  at  (7.2.  7.2.  7.2)  grid  units  with  0.5  kT 
total  energy  and  velocity  directed  in  z.  Compared  with  the  sphere  potential,  the 
electrons  are  cold. 


Table  4.1 


Grid  X  Y  2  Vx  Vy  Vz  Type  Energy  Comments 


1 

7.41 

7.21 

1  42 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

7.61 

7.21 

1  44 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

7  81 

7.21 

1.46 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

8.01 

7.21 

1.49 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

8.21 

7.21 

1.53 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

8.41 

7.21 

1.58 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

8.61 

7.21 

1.63 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

8.81 

7  21 

1.69 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

9.01 

7.21 

1.76 

0 

0 

1 

ELEC 

05 

Collected 

1 

9.21 

721 

1  84 

0 

0 

1 

ELEC 

0.5 

Collected 

1 

9.41 

7.21 

1.93 

0 

0 

1 

ELEC 

0.5 

Uncollected 

1 

9.61 

7.21 

2.03 

0 

0 

1 

ELEC 

0.5 

Uncollected 

1 

9.81 

7.21 

2  13 

0 

0 

1 

ELEC 

0.5 

Uncollected 

1 

10.01 

7  21 

2.25 

0 

0 

1 

ELEC 

05 

Uncollected 

1 

10.21 

7.21 

2.37 

0 

0 

1 

ELEC 

0.5 

Uncollected 

1 

10.41 

7  21 

2  51 

0 

0 

1 

ELEC 

0  5 

Uncollected 

1 

10.60 

7  21 

2.67 

0 

0 

1 

ELEC 

0.5 

Uncollected 

1 

10.80 

7.21 

2.83 

0 

0 

1 

ELEC 

05 

Uncollected 

1 

11.00 

7  21 

3.01 

0 

0 

1 

ELEC 

0.5 

Uncollected 

1 

11.20 

7.21 

3  22 

0 

0 

1 

ELEC 

0.5 

Uncollected 

be 

seen  in 

Fig.  4.4, 

the  first 

10  particles  were  collected  by  the 

object  while 

the  final  10  were  not.  This  corresponds  to  a  collection  radius  of  Rg  =  1.07  m. 
The  limiting  collection  distance  for  a  10-cm  sphere  at  11.000  volts  in  a  0.3  Gauss 
magnetic  field  given  by  Parker  and  Murphy  is  1.5  m. 
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4.3  Self-Consistent  Electron  Magnetosheath  Calculation 

The  procedure  to  compute  the  self-consistent  magnetosheath  described  in  Section  2 
was  used  for  the  10-cm  sphere  biased  to  11.000  volts.  The  calculation  of  the  particle 
trajectories  was  taken  up  following  the  calculation  of  the  space  charge  sheath  shown 
in  Fig.  4.3.  During  the  first  three  iterations  of  the  potential  trajectory  calculation. 
PWASON-CURREN,  the  number  of  left-right  pushes  in  z  was  limited  to  three.  This 
had  the  effect  of  limiting  the  lifetime  of  the  particles  in  the  sheath.  Consequently, 
the  shielding  due  to  these  particles  was  also  reduced,  creating  a  sheath  slightly 
broader  than  it  should  have  been.  This  in  turn  reduced  the  number  of  collectible 
orbits  and.  hence,  the  collected  current.  During  the  fourth  and  fifth  iterations,  the 
number  of  left-right  pushes  was  increased  to  five  and  six.  respectively.  In  Table  4.2. 
we  list  the  convergence  for  this  run. 

To  check  the  results,  we  developed  a  2  1/2-dimensional  computer  code  (2-dimensional 
coordinates  and  3-dimensional  velocities)  and  simulated  this  problem.  The  last  column 
in  Table  4.2  contains  the  results  of  this  code.  The  code  operates  in  cylindrical 

(r.  z.  0)  coordinates  on  a  quadrilateral  grid  consisting  of  90  spherical  shells  extending 

from  the  probe  surface  (0.1  m)  to  a  radius  of  7  m.  Each  shell  had  25  angular 

divisions  ranging  from  the  z-axis  to  the  r-axis.  The  z  =  0  plane  was  treated  as  a 

mirror  plane,  and  the  r  =  7m  sphere  was  held  fixed  at  zero  potential. 

Electrons  were  started  along  the  magnetic  field  (z)  direction  from  just  inside  the  7-m 
sphere  with  a  radial  spacing  of  2.3  cm.  and  with  current  and  velocity  set  to 
reproduce  the  one-sided  plasma  thermal  current  and  one-half  the  ambient  density. 
(Beyond  the  collection  radius,  the  remaining  half  of  the  ambient  density  comes  from 
electrons  reflected  from  the  z  =  0  mirror  plane.)  Energy  and  canonical  angular 
momentum  were  explicitly  conserved  at  each  particle  time  step.  Each  electron  was 
followed  until  it  either:  1.  hit  the  probe:  2.  escaped  the  grid;  or  3.  exceeded  the  limit 
of  how  many  times  it  was  allowed  to  bounce  off  the  mirror  plane.  At  each  time 
step,  the  particle’s  charge  was  bilinearly  shared  to  the  nearby  grid  points. 

_  u 

Ion  density  was  set  to  ne 

Iterations  consisted  of  alternately  solving  Poisson's  equation  and  tracking  particles  in 
the  resulting  potentials  to  update  the  charge  density.  This  process  was  continued 
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until  steady  state  solutions  were  achieved.  Potentials  were  constrained  to  be 
nonnegative;  and  in  most  runs,  charge  densities  were  constrained  to  be  nonpositive. 
(The  <.harge-density  constraint  affects  potentials  along  field  lines  within  the  collection 
radius  that  in  the  quasineutral  limit  satisfy 


M 

ne  =  0.5 


I*  had  negligible  effect  on  collected  currents  or  on  the  potentials  beyond  the  collection 

radius. 


Table  4.2 

Convergence  history  for  10-cm  sphere 


Iteration  Bounces  Current  (mA) 


2  1/2-D  Code  (mA) 


1  3  -4.93 

2  3 

3  3  -5.35 

4  5  -5.91 

5  6  -6.74 


4.2 

-6.7(10  bounces),  -11(50  bounces) 


As  can  be  seen  in  the  above  table,  the  current  in  the  POLAR  2.0  calculations  is  not 
♦i<r*v  converged.  Iteration  1  uses  potentials  computed  ignoring  magnetic  field  effects 
bu*  including  space  charge  effects.  It  compares  well  with  the  2  1/2-D  results.  The 
current  is  approaching  the  6.7  mA  result  predicted  by  the  2  1/2-D  self-consistent  code 
fo'  the  10-bounce  case.  The  sheath  current  for  Iteration  5  was  92  mA.  so  it  is 
evident  that  the  magnetic  field  is  dramatically  limiting  the  current  collection.  This  is 
expected  since  the  space  charge  sheath  without  a  magnetic  field  is  3.6  m  and  the 
maximum  collection  radius  given  by  canonical  angular  momentum  conservation,  the 
Parker  Murphy  limit,  is  1.5  m.  Hence,  this  calculation  is  in  the  strong  magnetic  field 
limited  regime 

Fig  4.5  shows  the  potentials  as  they  were  at  Iteration  5.  As  can  be  seen,  the 
contours  are  compressed  in  the  direction  perpendicular  to  the  magnetic  field,  which  is 
in  the  z  direction  This  is  due  to  particles  with  long  lifetimes.  The  sheath 
dimension  tries  to  adjust  to  fall  within  the  magnetic  limited  collection  radius. 
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5.  POLAR  2.0  SPEAR  I  CALCULATIONS 

This  section  contains  a  description  of  the  SPEAR  I  object  and  grid  structure  used  in 
the  POLAR  2.0  calculations.  Self-consistant  magnetosheath  calculations  are  presented 
for  two  SPEAR  I  voltage  configurations. 

5.1  SPEAR  I  Object  and  Grid  Definition 

The  POLAR  2.0  SPEAR  I  calculations  were  performed  in  the  object  and  computational 
grid  structure  shown  in  Fig.  5.1.  The  computational  grid  is  composed  of  a  total  of 
70.038  zones  in  four  grids.  The  three  inner  grids  provide  increased  resolution  for 
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the  rocket  body  (Grid  2).  the  plastic  boom  (Grid  3).  and  the  bushings  and  spheres 
(Grid  4).  Presently.  POLAR  2.0  does  not  allow  objects  to  cross  grid  boundaries  and 
does  not  allow  grids  to  share  an  edge.  This  restriction  results  in  gaps  in  the  SPEAR 
I  object  shown  in  the  figure.  These  gaps  are  small  and  should  not  cause  significant 
errors  in  the  calculations.  The  gap  between  the  bottom  of  the  boom  in  Grid  4  and 
the  top  of  the  boom  in  Grid  3  is  20  cm.  and  the  gap  in  the  rocket  body  between 

Grids  2  and  3  is  about  1  m. 

Grid  1 

Grid  1  is  the  outer  grid  and  is  designed  to  ensure  that  the  entire  space  charge  sheath 

will  be  contained  within  the  computational  space.  It  has  dimensions  34  by  30  by  30 

grid  units  and  a  mesh  size  of  0.5  m.  This  produces  a  computational  space  13.6  m 
by  12  m  by  12  m.  This  grid  contains  30.110  (=30.600-1792/64-29568/64) 
computational  *ones. 

Grid  2 

Grid  2.  the  lower  grid  in  Fig.  5.1.  contains  the  rocket  body.  Details  of  the  rocket 
body  are  shown  in  Fig.  5.2.  The  grid  is  8  by  8  by  28  Grid  2  units,  which  is  2  by  2 
by  6  Grid  1  units.  Grid  2  has  a  grid  ratio  of  four  with  respect  to  Grid  1,  producing 


SURFACE  cell  MATERIAL  COMPOSITION  AS  VIEWED  FROM  THE  POSITIVE  *  DIRECTION 
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a  mesh  size  of  0.1  m  and  a  computational  size  of  0.8  m  by  0.8  m  by  2.8  m.  the 
lower  portion  of  the  rocket  body  in  this  grid.  This  grid  contains  1.792  computational 
zones  2.6  m  long. 

Grid  3 

Grid  3  is  the  upper  grid  shown  in  Fig.  5.1.  which  contains  the  upper  rocket  body  and 
boom  structure.  This  grid  is  detailed  in  Fig.  5.3  and  the  upper  rocket  body  is 
detailed  on  Fig.  5.4.  The  grid  is  28  by  44  by  24  Grid  3  units,  which  is  7  by  11  by 
6  Grid  1  units.  This  grid  has  a  grid  ratio  of  four  with  respect  to  Grid  1.  producing 
a  mesh  length  of  0.1  m  and  a  size  of  2.8  m  by  4.4  m  by  2.4  m.  This  grid  contains 
Grid  4.  which  has  the  detailed  bushing  and  sphere  structures.  This  grid  contains 
28.344  (=29568-9792/8)  computational  zones. 

TMC  X  a  S  SLICE  OF  THE  CftlO  OF  THE  POL  Aft  PLOT  DATA  SET 
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SURFACE  CELL  MATERIAL  COMPOSITION  AS  VIEWED  FROM  THE  POSITIVE  X  DIRECTION 
FOR  X  VALUES  BETWEEN  1  AND  17 


MATERIAL  LECOtt 


i».  1«.  1«  14.  li.  #.  7.  t.  3.  1. 


Fig.  5.4 

Grid  4 

The  10-cm  spheres  and  resistive  grating  bushing  are  contained  in  Grid  4.  A  detailed 
view  of  this  grid  and  the  SPEAR  I  high  voltage  section  is  shown  in  Fig.  5  5.  This 
grid  is  34  by  24  by  12  Grid  4  units,  which  is  17  by  12  by  6  Grid  3  units.  This 
grid  contains  9,792  computational  zones. 
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5.2  SPEAR  I  Calculations 

POLAR  2.0  was  used  to  calculate  the  currents  to  the  SPEAR  I  rocket,  bushings,  and 
spheres  for  two  voltage  configurations.  The  first  configuration  was  a  34,000  volt- 
12.000  volt  combination  Down  the  bushings,  the  potentials  were  grated  smoothly  to 
zero  at  the  joining  point  of  the  two  supports  The  plasma  density  was  taken  to  be 
l.ell  m  with  a  temperature  of  0.1  eV.  The  second  calculation  was  similar  to  the 
first  case  except  that  the  voltages  were  scaled  down  to  8,900  volts  for  Sphere  1  and 
3,100  volts  for  Sphere  2. 

For  both  POLAR  runs,  the  starting  potential  configuration  was  computed  using  the 
analytic  form  for  the  charge  density,  p  =  p{$.  E).  As  discussed,  the  form  chosen 
gives  the  correct  density  for  symmetric  geometry  in  the  absence  of  a  magnetic  field. 
Using  this  analytic  approach  for  the  zeroth  stage  of  the  calculations  speeds  the 
convergence  of  the  later  iterations  by  starting  with  a  sheath  that  has  approximately 
the  correct  length  parallel  to  the  magnetic  field  and  slightly  enlarged  perpendicular  to 
the  magnetic  field  Following  this  first  calculation,  a  series  of  iterations  were 
performed  that  consisted  of  particle  tracking  to  determine  the  charge  density 
(CURREN)  and  solving  Poisson's  equation  to  determine  the  potentials  (PWASON) 
These  iterations  are  continued  until  the  self-consistent  sheath,  and  hence,  the  currents 
to  the  SPEAR  object  are  constant  Typically,  this  takes  5  to  10  iterations  depending 
on  the  specifics  of  the  calculation. 

Case  1  34.000-12.000 

To  improve  the  stability  of  the  calculation,  the  number  of  left-right  pushes  in  Z. 
bounces,  were  increased  as  the  calculation  progressed  This  had  the  effect  of  slowly 
reducing  the  width  of  the  sheath  in  the  perpendicular  direction  preventing  unstable 
nonconvex  sheath  boundaries  from  forming.  For  each  PWASON.  nine  space-charge 
iterations  were  performed  to  ensure  convergence  of  the  potentials  The  number  of 
bounces  along  with  the  total  current  to  the  object  in  each  grid  is  shown  below  in 
Table  5.1. 
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Table  5.1 

Convergence  History  for  34,000-12.000 


Iteration  Bounces  Grid 

Current 

1 

3 

2 

0.0000e+00  AMPS  Sun  3/60 

1 

3 

3 

0.0000e+00  AMPS 

1 

3 

4 

-2.5854e-02  AMPS 

2 

3 

2 

0.0000e+00  AMPS  Sun  3/60 

2 

3 

3 

0.0000e+00  AMPS 

2 

3 

4 

-2.2358e-02  AMPS 

3 

5 

2 

0.0000e+00  AMPS  Sun  3/60 

3 

5 

3 

0.0000e+00  AMPS 

3 

5 

4 

-3.2196e-02  AMPS 

5 

5 

2 

0.0000e+00  AMPS  Ridge  3200 

5 

5 

3 

0.0000e+00  AMPS 

5 

5 

4 

-3.6028e-02  AMPS 

6 

10 

2 

0.0000e+00  AMPS  Ridge  3200 

6 

10 

3 

0.0000e+00  AMPS 

6 

10 

4 

-4.7793e-02  AMPS 

7 

10 

2 

0.0000e+00  AMPS  Ridge  3200 

7 

10 

3 

0.0000e+00  AMPS 

7 

10 

4 

-4.6838e-02  AMPS 

8 

10 

2 

0.0000e+00  AMPS  Ridge  3200 

8 

10 

3 

-1.6428e-05  AMPS 

8 

10 

4 

-4.8471e-02  AMPS 

As  can  be  seen  from  the 

table,  the  current 

has  converged  and  remains  constant  for 

the  last  three 

iterations. 

The  last  column 

in  the  table  indicates  which  computer 

performed  the 

calculation. 

The  calculation 

was  initially  run  on  a  Sun  3/60,  and  then 

the  data  files  transferred  to  the  Ridge  3200  for  completion  of  the  calculations.  The 


breakdown  of  currents  to  the  object  is  given  below  in  Table  5.2. 
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Table  5.2 

Breakdown  of  Currents  for  34.000-12.000  Volt  Run 


Object  Current  (mA) 


Sphere  1 

-17.05 

Sphere  2 

-11.58 

Bushing  1 

-10.25 

Bushing  2 

-9.45 

Rocket  body 

0.00 

Plastic  booms 

-0.16 

Sphere  1  and  Bushing  1 

-27.30 

Sphere  2  and  Bushing  2 

-21.03 

The  total  current  is  -48.49  mA. 

The  total  current  of  48.49  mA  is  more  than  27  mA,  which  is  what  the  Parker- 
Murphy  limit  would  be  for  a  single  sphere  with  the  same  charge  and  an  average 
potential  of  23.000  volts.  This  is  due  to  breaking  of  the  symmetry  and  the  current 
to  the  bushings.  This  collected  current  is  much  smaller  than  the  sheath  current  for 
this  same  effective  sphere,  which  is  480  mA.  Fig.  5.6  shows  the  converged  potentials 
for  this  voltage  configuration. 


TMt  Z.  i  •  0000  4*0 TlWTWl  SUCf  FWOM  CYCi  41  Of  MTTHAK  77  Of  PO«SO* 


CMCT  Q*C  UWTT3 


Fig.  5.6  Converged  potentials  for  34.000-12.000  volt  run. 
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Ca«e  2  8,900-3,100 

Following  the  zeroth  calculation  using  the  analytic  form  for  the  charge  density,  six 
CURREN-PWASON  iterations  were  performed.  The  sequence  of  CURREN  runs  used 
the  following  number  of  bounces:  3.  3.  5,  10,  10.  For  each  PWASON.  nine 
space-charge  iterations  were  performed  to  ensure  that  the  potentials  were  converged. 

In  addition,  the  potentials  were  examined  after  each  iteration  to  determine  the  quality 
of  the  potential  solution.  The  convergence  of  the  total  current  collected  for  each  grid 
is  shown  in  Table  5.3  and  the  breakdown  of  the  currents  to  the  SPEAR  I  object  are 
shown  in  Table  5.4  below: 


Table  5.3 

Convergence  History  for  8,900-3.100 


Iteration 

Bounces 

Grid 

Current 

1 

3 

1 

0.0000e+00  AMPS 

1 

3 

2 

0.0000e+00  AMPS 

1 

3 

3 

-1.1168e-02  AMPS 

2 

3 

1 

0.0000e+00  AMPS 

2 

3 

2 

0.0000e+00  AMPS 

2 

3 

3 

-1.4087e-02  AMPS 

3 

5 

1 

0.0000e+00  AMPS 

3 

5 

2 

0.0000e+00  AMPS 

3 

5 

3 

-2.1766e-02  AMPS 

4 

10 

1 

0.0000e+00  AMPS 

4 

10 

2 

-2.1683e-06  AMPS 

4 

10 

3 

-3.0847e-02  AMPS 

5 

10 

1 

0.0000e+00  AMPS 

5 

10 

2 

0  0000e+00  AMPS 

5 

10 

3 

-2.7342e-02  AMPS 

6 

10 

1 

0.0000e+00  AMPS 

6 

10 

2 

0.0000e+00  AMPS 

6 

10 

3 

-3.2653e-02  AMPS 
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Table  5.4 

i 

i  Breakdown  of  Currents  for  8.900-3.100  Volt  Run 


Object  Current  (mA) 


Sphere  1 

-12.12 

Sphere  2 

-10.32 

Bushing  1 

-5.88 

Bushing  2 

-4.25 

Rocket  body 

0.00 

Plastic  booms 

-0.09 

Sphere  1  and  Bushing  1 

-18.00 

Sphere  2  and  Bushing  2 

-14.57 

The  total  current  is  -32.65. 

The  total  current  of  32.65  mA  is  more  than  14  mA.  which  is  what  the  Parker- 
Murphy  limit  would  be  for  a  single  sphere  with  the  same  charge  and  an  average 
potential  of  23.000  volts.  This  is  due  to  breaking  of  the  symmetry  and  the  current 
to  the  bushings.  This  collected  current  is  much  smaller  than  the  sheath  current  for 
this  same  effective  sphere,  which  is  156  mA.  Fig.  5.7  shows  the  converged  potentials 
for  this  voltage  configuration. 

THE  *-«•  0P00  POTEWTW.  SUC€  4*000  44  00  NTtRA*  ’0O»  PCMSSON 
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»  7, 

3i  3 1 

Ml, 
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f  4 
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00  41  M  144  1*3  Ml  M0  »T  Ml 
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Fig.  5.7  Converged  potentials  for  8.900-3.100  volt  run. 
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6.  COMPARISONS  POLAR  AND  NASCAP/LEO  SPEAR  I  CURRENT 
COLLECTION  CALCULATIONS 

Calculations  performed  using  both  the  POLAR  2.0  and  NASCAP/LEO  codes  showed 
basic  agreement  with  magnetic  limited  spherical  probe  theory  for  the  total  current. 

Ail  calculations  were  performed  for  ne  =  lO^m  3.  Te  =  O.leV.  and  b  =  0.4  Gauss. 

Table  6.1  shows  how  the  self-consistent  space  charge  effects  in  the  POLAR  2.0 
calculations  increase  the  currents  by  about  a  factor  of  two  over  the  NASCAP/LEO 
results.  Both  POLAR  2.0  and  NASCAP/LEO  predict  currents  dramatically  reduced  by 
the  magnetic  field:  for  example,  in  the  absence  of  magnetic  fields.  0.47  amperes  would 
be  collected  by  SPEAR  I  for  the  46,000  volt,  f  =  1/3  case.  The  initial  iteration  in 
the  POLAR  2.0  calculation  gave  total  currents  that  were  almost  identical  with  those 
from  NASCAP/LEO. 


Table  6.1 

Total  currents  in  amperes  collected  by  both  spheres  for  a  voltage  ratio  of  one-third. 
POLAR  2.0  includes  both  3-dimensional  geometry  and  self-consistent  space  charge. 


Voltage  PM  POLAR  NASCAP/LEO 


12.000  1.25e-2  3.26e-2  1.41e-2 

46.000  2.64e-2  4.85e-2  2.80e-2 


Table  6.2  shows  how  the  current  is  distributed  between  the  two  spheres  when  one 
sphere  has  about  three  times  the  potential  of  the  other  sphere.  For  the  lowest 
potential,  the  currents  are  divided  much  like  independent  probe  calculations  would 
predict.  In  the  NASCAP/LEO  calculations,  as  the  potentials  increase,  the  fraction 
going  to  the  higher  potential  sphere  increases.  We  attribute  this  to  an  overlap  of  the 
sheaths:  electrons  in  the  shared  region  are  preferentially  collected  by  the  higher 
potential  sphere.  However,  in  the  POLAR  2.0  calculations,  the  fraction  to  the  lower 
potential  sphere  initially  remains  near  the  independent  sphere  results.  As  the  magnetic 
field  compresses  the  sheath,  the  collection  by  the  two  spheres  becomes  more  even. 

This  difference  between  POLAR  2.0  and  NASCAP/LEO  has  not  been  tracked  down. 
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Both  particle  pushing  routines  have  been  tested;  the  difference  is  more  likely  due  to 
differences  in  the  computational  object  definition  or  small  differences  in  the  potential 
profiles.  More  work  is  needed  to  be  able  to  extrapolate  SPEAR  I  results  to  SPEAR 
II. 


Table  6.2 

Current  ratios  between  the  two  spheres  when  the  voltage  ratio  between  the  spheres  is 
one-third.  Values  listed  are  the  ratios  of  the  smaller  current  to  the  larger  current. 


Voltage  PM  POLAR  NASCAP/LEO 


12.000  0.59  0.81 

46.000  0.60  0.77 


0.533 

0.592 


7.  SUMMARY 

Further  calculations  are  planned  to  test  the  POLAR  2.0  code.  The  charge  module, 
electric  field  boundary  conditions,  and  wake  module  have  not  been  tested. 

Furthermore,  differences  were  found  between  the  NASCAP/LEO  and  the  POLAR  2.0 
codes  in  the  distribution  of  currents  to  the  SPEAR  I  rocket.  Because  of  the  large 
sheath  sizes  in  the  calculations,  these  differences  could  have  resulted  from  small 
differences  in  the  potentials  and  object  definitions.  This  needs  further  study. 

The  POLAR  2.0  support  code,  SHONTL,  is  not  completed.  Presently,  potential  plots 
and  tabular  data  can  be  displayed;  but  the  full  range  of  plots  has  not  been  fully 
implemented.  In  addition,  an  Iris  version  of  the  code  that  will  provide  convenient 
viewing  capability  is  being  prepared. 


